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We study charge and energy transfer in two-site molecular electronic junctions in which elec¬ 
tron transport is assisted by a vibrational mode. To understand the role of mode harmonic- 
ity/anharmonicity on transport behavior, we consider two limiting situations: (i) the mode is 
assumed harmonic, (ii) we truncate the mode spectrum to include only two levels, to represent 
an anharmonic molecular mode. Based on the models’ cumulant generating functions, we analyze 
the linear-response and nonlinear performance of these junctions and demonstrate that while the 
electrical and thermal conductances are sensitive to whether the mode is harmonic/anharmonic, 
the Seebeck coefficient, the thermoelectric figure-of-merit, and the thermoelectric efficiency beyond 
linear response, conceal this information. 


I. INTRODUCTION 

Molecular electronic junctions offer a rich playground for exploring basic and practical questions in quantum trans¬ 
port, such as the interplay between electronic and nuclear dynamics in nonequilibrium situations. Theoretical and 
computational efforts based on minimal model Hamiltonians are largely focused on the Anderson impurity dot model 
which consists a single molecular electronic orbital directly coupled to biased metal leads, as well as to a particular vi¬ 
brational mode [T]. Since the same molecular orbital is assumed to extend both contacts, the model allows simulations 
of transport characteristics in conjugated molecular junctions with delocalized electrons. 

In this work, we focus on a different class of molecular junctions as depicted in Fig. In such systems, two electronic 
levels are coupled via a weak tunneling element, but electrons may effectively hop between these electronic states when 
interacting with a vibrational mode. This could e.g., correspond to a torsional motion bringing orthogonal tt systems 
into an overlap as in the 2,2’- dimethylbiphenyl (DMBP) molecule recently examined in Refs. [2H1]- This model is 
related to the original Aviram-Ratner construction for a donor-acceptor molecular rectifier [S], thus we identify the 
two states here as D and A, see Fig. [l] and refer to the model as the “D-A junction”. More recently, this construction 
was employed for exploring vibrational heating and instability under a large bias voltage [SHE]- The system is also 
referred to as the “dimer molecular junction” [3], or an “open spin-boson model” [10] (where the spin here represents 
the D and A states, the bosons correspond to the molecular vibrational modes, and the system is open- coupled to 
metal leads). It was utilized to study charge transfer in donor-bridge-acceptor organic molecules |llj and organic 
molecular semiconductors as well as thermoelectric effects in quantum dot devices mm]. Recently, Erpenbeck 
et al. had provided a thorough computational study of transport characteristics with nondiagonal (or nonlocal) as 
well as diagonal (local) electron-vibration interactions |15| . Here, in contrast, we simplify the junction model and 
omit the contribution of direct tunneling between the D and A units. This simplification allows us to derive a closed 
(perturbative) expression for the cumulant generating function (CGF) of the model, which contains comprehensive 
information over transport characteristics. 

Measurements of charge current and electrical conductance in single molecules hand over detailed energetic and 
dynamical information m- Complementing electrical conductance measurements, the thermopower, a linear response 
quantity, also referred to as the Seebeck coefficient, is utilized as an independent tool for probing the energetics of 
molecular junctions, see for example Refs. mm]. Experimental efforts identified orbital hybridization, contact- 
molecule energy coupling and geometry, and whether the conductance is HOMO or LUMO dominated. More generally, 
the thermoelectric performance beyond linear response is of interest, with the two metal leads maintained at (largely) 
different temperatures and chemical potentials. 

What information can linear and nonlinear thermoelectric transport coefficients reveal on molecular junctions? 
Specifically, can they expose the underlying electron-phonon interactions and the characteristics of the vibrational 
modes participating in the process? Focusing on the challenge of efficient thermoelectric systems, how should we 
tune molecular parameters to improve heat to work conversion efficiency? These questions were examined in recent 
studies, a non-exhaustive list includes musnsHss]. 

We focus here on the effect of vibrational anharmonicity on thermoelectric transport within the D-A model. To 
explore this issue, two limiting variants of the basic construction are examined, as displayed in Fig. (a) The vibration 
is harmonic in the so-called “harmonic oscillator” (HO) model, (b) To learn about deviations from the harmonic 
picture, we truncate the vibrational spectrum to include only its two lowest levels, constructing the “anharmonic” 


2 


FIG. 1: Scheme of the D-A model considered in this work. A two-site (D, A) electronic junction is coupled to either (a) a 
harmonic molecular mode, or (b) an anharmonic mode, represented by a two-state system. The molecular mode may further 
relax its energy to a phononic thermal reservoir (not depicted here), maintained at temperature Tph. 



(AH) mode model. In a different context, the AH model could represent transport through molecular magnets, in 
which electron transfer is controlled by a spin impurity m- 

The complete information over transport behavior is contained in the respective cumulant generating functions, 
which we provide here for the HO and the AH models, valid under the approximation of weak electron-vibration 
interactions. From the CGFs, we derive expressions for charge and heat currents, and study the linear and nonlinear 
thermoelectric performance of the junctions. Focusing on the role of vibrational anharmonicity, we find that while it 
significantly influences the electrical and thermal conductances, nevertheless in the present model it does not affect 
heat-to-work conversion efficiency. 


II. MODEL 

We consider a two-site junction, where electron hopping between the D and A electronic states (creation operators 
cj; and cjj, respectively) is assisted by a vibrational mode. The total Hamiltonian is written as 

H ~ T T ^vih T ^el — vib 

+ eic\ci + erclcr- 

leL reR 

+ '^VlXclcd ++ '^Vr{clca + clcr). ( 1 ) 

IGL r^R 

The molecular electronic states of energies €d,a are hybridized with their adjacent metals, collection of noninteracting 
electrons, by hopping elements vi and v^. Here c| (cj) is a fermionic creation (annihilation) operator. The electronic 
Hamiltonian [Eq. Q, excluding Hyib + Hgi-yib] can be diagonalized and expressed in terms of new fermionic operators 
ai and a^. 


Hei ='^eialai+ '^eralar. ( 2 ) 

I r 

The molecular operators can be expanded in the new basis as, Cd = Ca = '^ylrO'r- The 7 /^^ coefficients 

satisfy 7 , = vi j (e/ - - Y.v ei-^+is ) The operators ci and can be expressed in terms of the new basis as 

Cl = with rjii' = Sill — ■ Similar expressions hold for the r set. 

Back to Eq. ([^, Hyib and Hei-yib represent the Hamiltonians of the molecular vibrational mode and its coupling 
to electrons, respectively. We assume an “off-diagonal” interaction with electron hopping between local sites assisted 
by the vibrational mode. Assuming a local harmonic mode we write 

Hyib — 

= 9[cW + cicd]{bl + bo), 


^el — vib 


(3) 
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with bo (6j) as the annihilation (creation) operator for the vibrational mode of frequency wq, g is the coupling 
parameter. The Hamiltonian Q then becomes 

Hho = Hei+uJoblbo + g ^ + /i.c.] (fej + 6o)- (4) 

leL,reR 

The second model considered here includes an anharmonic two-state mode. It is convenient to represent it with the 
Pauli matrices <Jx,y,zi and write the total Hamiltonian for the junction as 

Hah = Hei+ ^CT^ + g ^ [jijralar + h.c.]ax. (5) 

leL,reR 

The two models, Eqs. Q and ([^, describe electron-hole pair generation/annihilation by de-excitation/excitation 
of an “impurity” (vibrational mode). The left and right reservoirs defining H^i in Eq. (§ are characterized by a 
structured density of states since we had absorbed the D state in the L terminal, and similarly, the A level in the R 
metal. These electronic reservoirs are prepared in a thermodynamic state of temperature T^, = l/(fc_B/?,y) and chemical 
potential v = L,R^ set relative to the equilibrium chemical potential fip = 0. In our description we work with 
h = 1 and e = 1. Units are revived in simulations. 


III. TRANSPORT 

The complete information over steady state charge and energy transport properties of molecular junctions is 
delivered by the so-called cumulant generating function defined in terms of the characteristic function Z as 
5(Ae,Ap) = limt_>oo y ln2^(Ae, Ap) with Z{Xe,Xp) = (t))^ Here Ae,Ap are counting 

fields for energy and particles, respectively, defined for the right lead measurement, t is the final measurement time. 
The operators in this definition are Nh = the number operator for the total charge in the right lead, and 

similarly Hfi = epajap as the total energy in the same compartment. The superscript H identifies the Heisenberg 
picture, with operators evolving with respect to the full Hamiltonian. While closed results for the CGF can be derived 
for junctions of noninteracting particles [38], it is challenging to calculate this function analytically for models with 
interactions, see for example Ref. [39j . Our simplified D-A model is one of the very few many-body models which 
can be solved analytically. 

The CGF of the harmonic-mode junction Q can be derived using the nonequilibrium Green’s function (NEGF) 
technique [401141] assuming weak interaction between electrons and the particular vibration, employing the random 
phase approximation (RPA) |42L 143] . This scheme involves a summation over a particular set (infinite) of diagrams 
(ring type) in the perturbative series, taking into account all electron scattering processes which are facilitated by 
the absorption or emission of a single quantum wq. Physically, this summation collects not only sequential tunnelling 
electrons, but all coordinated multi-tunnelling processes, albeit with each electron interacting with the mode to the 
lowest order. The derivation of the GGF is nontrivial, and it is included in a separate communication |44| . Here we 
provide the final result 

Qho{K,Xp) = i(fcd - ku) - {ku + kdY - 4fc^)fc^. (6) 

The CGF of the AH model ([^ can be derived based on a counting-field dependent master equation approach [Hill], 

GAni^ej Ap) = —-{ku + kd) + 2 \/~ ^ ^u^d- C^) 

Both expressions are correct to second-order in the electron-vibration coupling g. It is remarkable to note on the 
similarity of these expressions which were derived from separate approaches. We use the short notation A = (Ap, Ae), 
where the counting fields are defined for right-lead measurements. It can be proved that our CGFs satisfy the 
fluctuation symmetry [45] 

G{Xe, Ap) = G{ — Xe — iAP, —Ap — — Pr^r)), (8) 

with A/3 = Pfi — /3i. This result is not trivial: Schemes involving truncation of interaction elements may leave out 
terms inconsistently with the fluctuation symmetry. Equations Q and Q are expressed in terms of an upward 
(excitation) and a downward (de-excitation) rates between vibrational states. The rates are additive in the two 
baths. 


( 9 ) 
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and obey the relation = k^[u}o —> — wo]- They are given by [71 Hi] 

Hf 

|^/L(e)(l - fnie + u;o))JL{e)JRie + ojo) e-*^.-^(^+-o)Ae^ 

/ OO 7 

^/fl(e)(l — /^(e + Wo)) Ji?,(e) JL(e + wq) . (10) 


The rates k^ and k^ are evaluated from these expressions at A = 0; f^(e) = [exp(/3^(e — fi^)) + 1]“^ is the Fermi- 
Dirac distribution function of the v = L,R lead. The properties of the molecular junction are embedded within 
the spectral density functions, peaked around the molecular electronic energies ed,a with the broadening ri.(e) = 


Jhii) 

Tfl(e) 


9 
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rL(e) 

(e-e,)2 + (rZ(e)/2)2’ 

rfl(e) 

(e-ea)2 + (r«(e)/2)2' 


( 11 ) 


These expressions are reached through the diagonalization procedure of the electronic Hamiltonian while ignoring the 
real principal value term, responsible for a small energy shift of e^.a [7|. In what follows we take T^ as a constant 
independent of energy and assume broad bands with a large cutoff the largest energy scale in the problem. 

We obtain currents and high order cumulants by taking derivatives of the CGF with respect to the counting fields. 
The particle (Ip) and energy {1^} current are given by 


(4)^ 


t 


( 4 ) = 


M 

t 


dg{Xe,Xp) 

d{iXp) 


Ae 


— Ap—0 


dg{Xe,Xp) 

d{iXe) 


Ae —Ap — 0 


( 12 ) 


After some manipulations we reach the compact form for the harmonic (—) and anharmonic (+) models. 


^iAh/ho^ = 2 


uR — vLt^R— vL kL — vRt^L— 

kd ± ku 

kd [9(ae)*«lA.=o] 


»[%AjfcdU.=o] 


kd ± kn 


(13) 


The rates are given by Eq. (10) with A = 0. It is notable that the only difference between the HO and AH models is 
the sign in the denominator. Note that we did not simplify the expression for the energy current above; 

the derivatives return energy transfer rates analogous to Eq. (10), only with an additional energy variable in the 
integrand. While figures below only display quantities related to charge and energy currents, it is useful to emphasize 
that the CGF contains information on fluctuations of these currents. For example, the zero-frequency noise for charge 
current is given from 


(Sp)^ 


{{N^}) d^g{Xe,Xp) 


d{iXp)‘^ 


Ae — A ry — 0 


(14) 


where {{N"^)) = {N'^) — {N)'^ is the second cumulant. 

Our derivation is based on the diagonal representation of the electronic Hamiltonian, thus the occupations of 
the molecular electronic states D and A follow the Fermi function by construction. In the weak coupling limit 
employed here, the back-action of the vibrational degrees of freedom on the electronic distribution is not included. 
While in other models [IS] this back-action may be significant, here we argue that its role is rather small: Recent 
numerically-exact path integral simulations [8] testify that this type of quantum master equation very well performs 
at weak-intermediate electron-vibration coupling, justifying our scheme. Note that in path integral simulations [S] 
the states D and A were absorbed into the metal leads as well, yet the electronic distribution was allowed to evolve 
in time, naturally incorporating the back-effect of vibrations on the electronic distribution in the steady-state limit. 


IV. RESULTS 

We are interested in identifying signatures of mode harmonicity in transport characteristics. We set the right 
contact as hot, Tr > Tr, and write the electronic heat current extracted from the hot bath by (Ih) = (4) — k-R{Ip)- 
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FIG. 2: Linear response behavior of the D-A junction as a function of molecule-metal hybridization with a harmonic mode 
(full) and an anharmonic two-state system mode (dashed), (a) Normalized electrical conductance G/Go with Go = ^ jh, the 
quantum of conductance per channel per spin, (b) Seebeck coefficient S. (c) Electronic thermal conductance S, and (d) the 
figure of merit ZT. Parameters are eo = 0.01, ujq = 0.02, g — 0.01 in eV, room temperature T = 300 K. We assumed flat bands 
with a constant density of states. 






The bias is applied such that > /ifl, thus the macroscopic efficiency of a thermoelectric device, converting heat to 
work, is given by 


"‘IW- ih) ^ 

The device is operating as a thermoelectric engine when both charge and energy current flow from the hot (right) 
bath to the cold one. Note that according to our conventions the currents are positive when flowing from the right 
contact to the left. 


A. Linear response coefficients 


In linear response, i.e., close to equilibrium, the charge current (Ip) and heat current (Ih) as obtained from Eq. (13) 


can be expanded to lowest order in the bias voltage Afj, = — hl = eV and temperature difference AT = Tr — Tr. 

To re-introduce physical dimensions, we multiply the charge current by ejh and the heat current by 1/h. The resulting 
expansions are cumbersome thus we write them formally in terms of the coefficients Oij, (i,j = h,p), 


6 6 

{Ip} — j^cLppApT j^dp fik^AT^ 

ih) = ^ah,pAp+ ^ah^hkBAT. 


(16) 


Recalling that {Ip) = GV + GSAT and {Ih) = GTW -b (E -|- G'S'n)AT with 11 as Peltier coefficient [371 EH]: we 
identify the electrical conductance G = \ap,p, the thermopower S = the electron contribution to the thermal 

^ and the (dimensionless) thermoelectric figure of merit ZT = ^^§-T which 


conductance E = ^ 


ah,h 


determine the linear response thermoelectric efficiency. We obtain these coefficients numerically, by simulating Eq. 
(13) under small biases. 

^ display the behavior of G 


Figures |2]|^ display the behavior of G, S, E and ZT at room temperature T = 300 K for the harmonic and 
anharmonic-mode junctions. In numerical simulations below the phononic contribution to the thermal conductance 
is ignored, assuming it to be small compared to its electronic counterpart. A quantitative analysis of the contribution 
of the phononic thermal conductance is included in the Discussion Section. In addition, for simplicity, the junction 
is made spatially symmetric with T = Tl,r and eo = ed,a- The currents are given by Eq. ( [l^ , and we make the 
following observations: (i) The harmonic-mode model supports higher currents relative to the two-state case, but 
at low temperatures, wo/T))l, when the excitation rate is negligible relative to the relaxation rate, the two models 
provide the same results, (ii) Since the expressions for the currents in the HO and the AH models are proportional 
to each other, the resulting thermopower and figure of merit are identical. 
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FIG. 3: Linear response behavior of the donor-acceptor junction as a function of gate voltage, (a) 
Seebeck coefficient, (c) electronic thermal conductance, and (d) figure of merit ZT. Parameters 


r = 0.05 eV. 


Electrical conductance, (b) 
are the same as Fig. for 



Eq [eV] 


eo[eV] 


FIG. 4: Linear response behavior of the donor-acceptor junction as a function of vibrational frequency cjo for eq = 0.2 eV, 
F = 0.2 eV, g — 0.01 eV, and T=300 K. (a) Electrical conductance, (b) Seebeck coefficient, (c) electronic thermal conductance, 
and (d) figure of merit ZT. 



Figure [^displays transport coefficients as a function of metal-molecule hybridization assuming a resonance situation 
ksT > eg. The conductances show a turnover behavior in accord with Eq. 0’ growing with F for small values 
F < eo, then falling down approximately as G, S oc F”^. The figure of merit shows a monotonic behavior, increasing 
when levels’ broadening becomes small r((r as we approach the so called “tight coupling” limit in which charge and 
heat currents are (optimally) proportional to each other. ZT can be significantly enhanced by tuning the molecule to 
an off-resonance situation, eg > ksTjF, see Fig. We find that the electrical and thermal conductances strongly fall 
off with cq) but the Seebeck coefficient displays a non-monotonic structure, with a maximum showing up off-resonance 
[49j . resulting in a similar enhancement of ZT around eg = 0.2. It can be proved that the conductances are even 
functions in gate voltage, G(eo) = G(—cq), E(eo) = E(—eo) while S'(eo) = —S(—eo), resulting in an even symmetry 
for ZT with gate voltage. 

In Fig. we show transport coefficients as a function of the vibrational frequency. Parameters correspond to a 
resonant situation eg/V = 1. Both G and E decay exponentially with ojg when ojg > ksT. However, the figure of 
merit only modestly increases with ojg in the analyzed range due to the enhancement of S in this region. The values 
reported for ZT in figure]^ can be increased by weakening the metal-molecule coupling energy F. 

The maximal efficiency, r]max = and the efficiency at maximum power, ri{Pmax) = ^ zt '+2 are 

shown in Fig. as a function of eo and F for a fixed molecular frequency ojg = 0.02 eV and temperature T = 300 
K. Here, ?7c = 1 — Tcoid/Thot corresponds to the Carnot efficiency. By tuning the gate voltage and the molecule-lead 


























7 


FIG. 5: Contour plot of linear response efficiencies as a function of hybridization F and electronic energies (or gate voltage) cq. 
(a) Maximum efficiency, (b) Efficiency at maximum power. Parameters are cuo = 0.02 eV and T= 300 K. 

(a) 


ri /ri^ 
'max 'C 


(b) n(P,J/Tlc 


> 

u 



0.1 0.2 0.3 



0.1 0.2 0.3 


EoIeV] 


Eq [eV] 


FIG. 6: Transport beyond linear response, (a) current voltage cha racteristics for the harmonic (full) and anharmonic (dashed) 
mode models, (b) Heat to work conversion efficiency [Eq. (15l]. Parameters are cjo = 0.02, eo = 0.2, g=0.01, P = 0.1, 
Pph = 0.002 in units of eV, and Tl= 300 K, Ti{= 800 K and rpji=300 K. 



hybridization we approach the bounds rjmax/vc 1) v{Pmax)/VC 1/2 [ll]. Particularly, for V ^ 0.01 eV we receive 

Vmax/ric = 0.8 at the energy eo = 0.15. 


B. Nonlinear performance 


Nonlinear thermoelectric phenomena are anticipated to enhance thermoelectric response |50| . Elastic scattering 
theories of nonlinear thermoelectric transport have been developed e.g. in Refs. accounting for many-body 

effects in a phenomenological manner. Only few studies had considered this problem with explicit electron-phonon 
interactions, based on the Anderson-Holstein model [55j or Fermi-Golden rule expressions [la¬ 
in Fig. we simulate the current-voltage characteristics and the resulting efficiency of the D-A junction beyond 
linear response, by directly applying Eq. (13). As discussed in previous investigations EHS], the molecular junction 
may break down far from equilibrium due to the development of “vibrational instability”. This over-heating effect 
occurs when (electron-induced) vibrational excitation rates exceed relaxation rates. To cure this physical problem, 
we allow the particular vibrational mode of frequency uq to relax its excess energy to a secondary phonon bath of 
temperature Tph- This can be done rigorously at the level of the quantum master equations and within the NEGF 
technique [Hill] to yield the rates ka = k^^’^+ +Tph{ujQ)[nph{ciQ) + 1], = k^~^^ + k^^^ + Tph{ijJo)nph{ujQ), 

with UphiuJo) = [e‘^o/kBTph _ ^j-i ^ damping term rp;i(a;o). Interestingly, we confirmed (not shown) that this 
additional energy relaxation process does not modify the thermoelectric efficiency displayed in Fig. ib). 

















V. DISCUSSION AND PROSPECT 


We focused on two-site electronic junctions in which electron transfer between sites is assisted by a particular 
mode, harmonic or anharmonic (two-state system). The complete information over steady state transport behavior 
is catered by the cumulant generating function, which we provide here for the HO and the AH mode models, valid 
under the approximation of weak electron-vibration interaction. We explored linear-response properties, the electrical 
and thermal conductances G and E, as well as the Seebeck coefficient S', the thermoelectric figure of merit ZT^ and 
the resulting efficiency. We further examined current-voltage behavior and the heat-to-work conversion efficiency 
far-from-equilibrium. We found that G and E (more generafiy, the charge and energy currents) are sensitive to the 
properties of the mode, while S and ZT are insensitive to whether we work with a harmonic mode or a truncated 
two-state model. Several comments are now in place: 

(i) Genuine anharmonicity. We examined the role of mode anharmonicity by devising a two-state impurity model. 
It should be emphasized that in the context of molecules, the two-state impurity does not well represent vibrational 
anharmonicity at high temperatures, as many states should then contribute. Furthermore, it misses an explicit 
parameter tuning the potential anharmonicity. However, the AH model allows for a first indication on how deviations 
from harmonicity reflect in transport behavior. The HO and AH models have similar CGFs, yielding currents which 
are proportional, thus an identical thermoelectric efficiency. It can be readily shown that an n-state truncated HO 
provides a figure of merit identical to the infinite-level HO model, but it is interesting to perform more realistic 
calculations and consider e.g. a morse potential to represent a physical anharmonic molecular vibration. In this case, 
an analytical form for the CGF is missing, but one could still derive the charge current directly from a quantum 
master equation formalism, to obtain the performance of the system. We expect that with a genuine anharmonic 
potential, S and ZT would show deviations from the harmonic limit, as different pathways for transport open up. 
Overall, we believe that our results here indicate on the minor role played by mode anharmonicity in determining 
heat-to-work conversion efficiency. 

(ii) Direct tunneling. Our analysis was performed while neglecting direct electron tunneling between the D and 
A sites. This effect could be approximately re-instituted by assuming that coherent transport proceeds in parallel 
to phonon-assisted conduction, accounting for the coherent contribution using a Landauer expression, see e.g.. Refs. 
[3111]. Indeed, path integral simulations indicted that in the D-A model, coherent and the incoherent contributions 
are approximately additive [5|- 

(hi) Strong electron-phonon interaction. The GGFs (§-0 are exact to all orders in the metal-molecule hybridization 
but perturbative (to the lowest nontrivial order) in the electron phonon coupling g. This is evident from the structure 
of the rate constants in Eq. (10 1 , as electron transfer is facilitated by the absorption/emission of a single quantum 
loq. In numerical simulations we typically employed g = 0.01 eV and wq = 0.02 eV. This value for g may seem large 
given the perturbative nature of our t reat ment requiring g/uJo{{l. However, since in the present weak-coupling limit 
the current simply scales as 5 ^, Eq. (13), our simulations in this work are representative, and can be immediately 


translated for other values for g. It is of interest to generalize our results and study the performance of the junction 
with strong electron-phonon interaction, e.g. by using a polaronic transformation |5bHfi0j . 

(iv) Phononic thermal conductance. We studied here electron transfer through molecular junctions, but did not 
discuss phonon transport characteristics across the junction, mediated by molecular vibrational modes. Gonsideration 
of the phononic thermal conductance Uph is particularly important for a reliable estimate of ZT, as the thermal 
conductance E should include contributions from both electrons and phonons. We now estimate Kph. The quantum of 
thermal conductance, an upper bound for ballistic conduction, is given by kq = TrkgT/GH |61j . At room temperature, 
this provides kq = 0.28 nW/K, which exceeds the electronic thermal conductance obtained in our simulations, to 
dominate the total thermal conductance and predict (significantly) lower values for the figure of merit. However, one 
should recognize that at high temperatures the ballistic bound for phonon thermal conductance is far from being 
saturated as was recently demonstrated in Ref. [62) . In particular, the phononic thermal conductance of a two-level 
junction was evaluated exactly in Ref. [53 j it significantly falls below the harmonic bound |52|- 

Eor a concrete estimate, we adopt the perturbative (weak mode-thermal bath) expression for the phononic current 
through a two-state junction developed in Ref. [64] and further examined in Ref. [65] - 


■AH 

Jph ~ 9ph^0 


[nR^uJo) - ULjuJo)] 

[2 -I- 2n/{(wo) -f 2ni(wo)] 


(17) 


Here, gph is the interaction strength of the local vibrational mode to the phononic environments at the two terminals 
(assuming identical interaction strengths). The baths are characterized by their Bose-Einstein distribution functions 
ni,(ijj). In the case of a local harmonic mode, Eq. 0 holds, only missing its denominator. Using wq = 0.02 eV, 
we receive an estimate for the phononic thermal conductance ~ 3.5 X gph nW/K. Thus, 

as long as the mode-bath coupling Ppp^ is taken as weak, for example, gph < 5 meV for the data of Fig. the 
electronic contribution to the thermal conductance dominates the total thermal conductance and our simulations are 
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intact. Similar considerations hold for harmonic-mode junctions. Proposals to reduce the coherent phononic thermal 
conductance by quantum interference effects [66] and through-space designs m could be further considered. 

(v) High order cumulants. The cumulant generating functions, Eqs. ([^, and Q, contain significant information. 
For example, one could examine the (zero-frequency) current noise, to find out to what extent it can reveal microscopic 
molecular information. 

(vi) Methodology development. The cumulant generating function of the HO model was derived from air NEGF 
approach |44| . The corresponding function for the AH model was reached from a master equation calculation IZlIll]. 
Both treatments are perturbative to second order in the electron-vibration interaction; we take into account all electron 
scattering processes that are facilitated by the absorption or emission of a single quantum ujq. It is yet surprising to 
note on the direct correspondence between NEGF and master equation results, as derivations proceeded on completely 
different lines. In particular, the NEGF approach was done at the level of the RPA approximation to guarantee the 
validity of the fluctuation theorem. The master equation approach has been employed before to study currents (first 
cumulants) in the HO model (3, showing exact agreement with NEGF expressions presented here. This agreement, 
as well as supporting path integral simulations [8], indicate on the accuracy and consistency of the master equation 
in the present model. Given its simplicity and transparency, it is of interest to extend this method and examine 
higher order processes in perturbation theory, to gain further insight on the role of electron-vibration interaction in 
molecular conduction. 

(vii) Efficiency fluctuations. We focused here only on averaged-macroscopic quantities. However, in small systems 
fluctuations in input heat and output power are significant, resulting in “second law violations” as predicted from the 
fluctuation theorem |45| . to e.g. grant efficiencies exceeding the thermodynamic bound. To analyze the distribution 
of efficiency, the concept of “stochastic efficiency” has been recently coined and examined [68l [69]. In a separate 
contribution [33] we extend the present analysis and describe the characteristics of the stochastic efficiency in our 
model, particularly, we explore signatures of mode anharmonicity in the statistics of efficiency. 

(viii) Molecular calculations. It is of interest to employ our expressions and examine heat-to-work conversion 
efficiency in realistic molecular junctions. Our results demonstrate that the conversion efficiency can be improved by 
working in the off-resonant limit, r/eo((l, as well, when tuning cq through a gate voltage to Co/ksT ~ 5. In such 
situations, the figure of merit ZT can be made large since one can tune S to large values (though the conductances 
are small). This suggests that in the D-A class of molecules one should focus on enhancing the thermopower as a 
promising mean for making an overall improvement in efficiency. 

In our ongoing work we are pursuing some of these topics. The derivation of the GGFs employed here and the 
behavior of efficiency fluctuations in linear response, and beyond that, are detailed in Ref. [33] ■ In Ref. m we 
derive thermoelectric transport coefficients for the dissipative D-A model, beyond linear response, and describe the 
operation of thermoelectric diodes and transistors. 
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